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Abstract 

A research project is underway at the NASA Glenn Research Center (GRC) to produce computer 
software that can accurately predict ice growth under any meteorological conditions for any aircraft 
surface. This report will present results from two different computer programs. The first program, 
LEWICE version 3.2.2, has been reported on previously. The second program is GlennICE version 0.1. 
An extensive comparison of the results in a quantifiable manner against the database of ice shapes that 
have been generated in the GRC Icing Research Tunnel (IRT) has also been performed, including 
additional data taken to extend the database in the Super-cooled Large Drop (SLD) regime. This paper 
will show the differences in ice shape between LEWICE 3.2.2, GlennICE and experimental data. This 
report will also provide a description of both programs. Comparisons are then made to recent additions to 
the SLD database and selected previous cases. Quantitative comparisons are shown for horn height, horn 
angle, icing limit, area and leading edge thickness. The results show that the predicted results for both 
programs are within the accuracy limits of the experimental data for the majority of cases. 

Introduction 

The Icing Branch at GRC has produced several computer codes (refs. 1 to 35) over the last twenty 
years for performing icing simulation. While some of these tools have been collaborative projects, most 
have been developed primarily by one person with some assistance by others. The state of computing has 
also changed dramatically in that time period. Early on, design decisions were made in these tools that 
took into consideration the limited capabilities of early computers. Additionally, the personnel tasked to 
design and write these computer codes were and remain icing experts, not software experts. As these 
codes have grown in complexity and have been accepted by users as production icing tools, there has 
arisen a need for the developers to adhere to standard software practices used to develop commercial 
software. The goals of these practices, as they relate to the GlennICE effort, include: creating a team- 
based development process rather than the current individual-based process; building modular code 
elements that can be re-used as needed by several processes; designing or redesigning code elements with 
future capabilities in mind; and developing efficient testing methodologies for testing and validation of 
the code. This report will, in part, provide an update on the status of that process and describe the current 
capabilities of the GlennICE software (refs. 36 and 37). 

This report will address the validation of this software against a recent set of ice shape data in the 
Super-cooled Large Droplet (SLD) regime (ref. 38). This validation effort mirrors a similar effort 
undertaken previously for previous validations of LEWICE (refs. 39 and 40). Those reports quantified the 
ice accretion prediction capabilities of the LEWICE software. Several ice geometry features were 
proposed for comparing ice shapes in a quantitative manner. The resulting analysis showed that LEWICE 
compared well to the available experimental data. 

The effects of super-cooled large droplets in icing have been researched extensively since the 1994 
ATR-72 crash in Roselawn, Indiana (ref. 41). It has been speculated that the accident occurred due to the 
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accumulation of ice aft of the deicing boots. This ice accumulation formed aft of the boots due to 
impingement of drops greater than 40 pm in diameter. Since then, several experimental efforts have been 
made to document SLD ice shapes and to investigate the underlying physics (refs. 42 to 49). While this 
report will provide comparisons to standard icing conditions, the emphasis was placed on the newer data, 
which is predominately SLD. 

The report is divided into four sections. The first section will provide a description of the LEWICE 
(ref. 34) software, the Smagglce (ref. 35) software and the Navier-Stokes flow solver WIND (ref. 50) that 
was used for this effort. The second section will provide a description of the GlennICE model with 
emphasis on the differences between the two development projects. The third section will describe the 
experimental data and the parameters used for quantifying the comparisons. The last section will provide 
validation results along with a statistical comparison of those parameters with the available experimental 
data. 


Nomenclature 

B ice height (m) 

c chord (m) 

c n series term in analytic temperature solution (m _1 ) 

c p specific heat (J/kgK) 

H heat transfer coefficient (W/m 2 K) 

h film thickness (m) 

k thermal conductivity (W/mK) 

K l LEWICE splashing parameter (dimensionless) 

L characteristic length (m) 

LWC liquid water content (g/m 5 ) 

m water mass flux (kg/m 2 s) 

n unit normal (dimensionless) 

P pressure (N/m 2 ) 

Q flow rate (m 3 /s) 

q heat flux (W/m 2 ) 

s surface wrap distance (m) 

T temperature (K) 

t time (s) 

V velocity (m/s) 

W impact velocity (m/s) 

w f wetness factor (dimensionless) 

v chordwise direction (m) 

Y airfoil metal thickness (m) 

y direction normal to the surface (m) 


Dimensionless numbers: 


c / 


skin friction coefficient = 


2T y=0 

PairVl 


Nu 


Nusselt number = 


HL 

k 
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Oh 

Ohnesorge number = . ^ 

yjpad 

Pr 

c p 

Prandtl number = — — 

k 

Re 

pVL 

Reynolds number = 

P 

St 

o , H 

Stanton number = 

Pair c pVo o 

Greek letters: 

a 

thermal diffusivity (m2/s) 

A 

latent heat (J/kg) 

P 

viscosity (kg/m*s) 

V 

kinematic viscosity of air (nr 

P 

density ( kg/m 3 ) 

Subscripts: 

a 

air 

acc 

accumulation 

cond 

conduction 

d 

drop 

e 

edge of boundary layer 

eff 

effective 

evap 

evaporation 

f 

freezing 

film 

film 

i 

ice 

imp 

impinging 

k 

roughness 

m 

mass transfer 

mp 

melting point 

o 

total 

rb 

runback 

rec 

recovery 

s 

surface 

splash 

splash 

t 

turbulent 

V 

vapor 

w 

water 

X 

x-dependent 

y 

y-dependent 

00 

free-stream property 
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Existing Computational Tools 

LEWICE 


The computer program, LEWICE, embodies an analytical ice accretion model that evaluates the 
thermodynamics of the freezing process that occurs when super-cooled droplets impinge on a body. The 
atmospheric parameters of temperature, pressure, and velocity, and the meteorological parameters of 
liquid water content (LWC), droplet diameter, and relative humidity are specified. This information, 
along with the body geometry, is used to determine the shape of the ice accretion through an evaluation of 
the mass and energy balance associated with freezing of the incoming liquid water. The surface of the 
clean (un-iced) geometry is defined by segments joining a set of discrete body coordinates. The software 
consists of four major modules. They are: (1) the flow field calculation, (2) the particle trajectory and 
impingement calculation, (3) the thermodynamic and ice growth calculation, and (4) the modification of 
the current geometry by addition of the ice growth. 

LEWICE applies a time-stepping procedure to “grow” the ice accretion. Initially, the flow field and 
droplet impingement characteristics are determined for the clean geometry. The ice growth rate on each 
segment defining the surface is then determined by applying the thermodynamic model. When a time 
increment is specified, this growth rate can be transformed into an ice thickness and the body coordinates 
are adjusted to account for the accreted ice. This procedure is repeated, beginning with the calculation of 
the flow field about the iced geometry, then continued until the desired icing time has been reached. The 
results shown in this report are from version 3.2.2 of LEWICE (ref. 34). 

WIND 

WIND-US (ref. 51) is an unstructured, multi-zone, compressible flow solver with flexible chemistry 
and turbulence models. Zonal interfaces may be abutting or overlapping, allowing the flexibility to treat 
complex systems moving relative to one another. WIND-US is a computational platform that may be used 
to numerically solve various sets of equations governing physical phenomena. Currently, the software 
supports the solution of the three-dimensional Euler and Navier-Stokes equations of fluid mechanics, 
along with supporting equation sets governing turbulent and chemically reacting flows. 

WIND-US is a product of the NPARC Alliance, a partnership between GRC and the Arnold 
Engineering Development Center (AEDC) dedicated to the establishment of a national, applications- 
oriented flow simulation capability. The Boeing Company has also been closely associated with the 
Alliance since its inception, and represents the interests of the NPARC User's Association. 

Smagglce 

The Surface Modeling and Grid Generation for Iced Airfoils (Smagglce) (ref. 35) toolkit is a suite of 
interactive tools developed at GRC which can simplify and improve an icing effects study. This suite of 
tools is used to prepare two-dimensional cross sections of iced airfoils for computational fluid dynamics 
(CFD) analysis. Smagglce is designed to help researchers and engineers study the effects of ice accretion 
on airfoil performance, which is difficult to do with other software packages because of the complexity of 
ice shapes (ref. 52). CFD tools are used primarily for certification studies and analysis, to evaluate safety 
of flight with unprotected surfaces and to determine the need for ice protection. In this study, Smagglce 
was used to generate grids for each airfoil and for interim ice shapes in the multi-time step accretion 
process. 


GlennICE 


The GlennICE project was initiated in 1999 as a means to develop a consistent set of standards and 
practices for creation, use and maintenance of NASA administered icing computational tools. One goal of 
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this process is to provide a mechanism for a more team-based approach to software development. Current 
development of icing codes at NASA has been more separated and individual-based. The plan has been to 
establish a common computational framework for integration of NASA icing analytical tools, and to 
create a set of requirements both for the development of these tools and for their integration into the 
common framework. The features of the resulting system would include providing a common user 
interface for all tools and allowing user selection of analysis modules. A modular design would allow 
more streamlined substitution of new or improved software elements, and the defined process would 
include documented standards for rapid integration of new software. This process also requires a version 
control system designed to allow both control of approved software and flexibility for rapid prototyping, 
as well as documentation that allows rapid identification of software design and use. The first step of this 
process has been to integrate existing icing computational tools and to redesign them if necessary so that 
they fit into the design. Where applicable, new or modified routines will be noted in this report. The 
following discussion will present each of the main subdivisions of the software and describe the 
modifications made. 


Grid Generation/Surface Geometry Modeling 

The methodology for generating surface points (panels if the potential flow solver is used) has been 
completely replaced. The current methodology relies on using the Smagglce software described earlier in 
order to generate the surface points and grids (if needed) for 2D geometries. The philosophy in LEWICE 
was that the user should not have much control over the point distribution in order to provide uniformity 
of results when run by users with different experience levels. Using Smagglce gives the user much more 
control over the point distribution and grids than could be done in LEWICE. Although Smagglce is a 
GUI-driven interactive software, the program can be run with scripts that do not require user interaction. 
This allows both flexibility in point distribution and automation for multi-time step runs. Prior to the 
release of GlennICE, it will be necessary to test this feature in order to ensure that users with little or no 
experience can produce reliable results. The LEWICE routines for automatic point distribution may be 
added to GlennICE in the future based on feedback from the users. In 3D, grids are input in the same 
format as the LEWICE3D (ref. 32) code. The grid formats supported by LEWICE3D can also be used for 
2D simulations. This feature allows the user additional flexibility in selecting a grid generator and flow 
solver since LEWICE3D supports a wide variety of formats. Automatic time stepping can only be used in 
this instance if the grid generator selected by the user can automatically produce accurate grids for the 
new geometry. 


Flow Solvers 

The potential flow solver in LEWICE was converted to Fortran90 so that it could be used in 
GlennICE. In addition, GlennICE can support both multi-block structured grids (max. 50 blocks) and 
unstructured grids. Additional integration with Navier-Stokes solvers was performed to use more of the 
flow solution results in the heat transfer routines. In LEWICE, the user could have heat transfer 
coefficients (or heat fluxes) read from a file. It was presumed that the heat transfer coefficients were 
calculated in the Navier-Stokes solver. Additionally, LEWICE would determine velocities at the edge of 
the boundary layer through use of Bernoulli’s equation from the surface pressures. While this method 
gave heat transfer results very similar to potential flow, it did not use the Navier-Stokes results to their 
fullest extent. In GlennICE, the heat transfer coefficient process was modularized so that the LEWICE 
heat transfer correlation could be used independently of the integral boundary layer. The LEWICE 
turbulent heat transfer correlation is taken from Schlichting (ref. 53) as 
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— c^Re x Pr 

Nu x = -2 — t — — where St k = 0.52 Re/' 45 Pr 0 ' 8 

Pr r + ^ c /% 

according to Owen and Thompson (ref. 54). This equation relies on the skin friction coefficient (available 
from the Navier-Stokes solver), the local Reynolds number and the roughness Reynolds number. Of these 
three values, only the roughness Reynolds number is not output by Navier-Stokes codes. The roughness 
Reynolds number can be calculated from the roughness height and the velocity at the roughness height 
estimated from the Navier-Stokes solution. In order to complete the analysis, the laminar heat transfer 
coefficient and transition location are needed as well. The laminar heat transfer coefficient is simply 
related to the local Reynolds number by 


1 V V 

Nu = — Pr /3 Re/ 2 
3 

Transition is determined when the roughness Reynolds number is greater than 600. This differs 
slightly from the value LEWICE uses. LEWICE varies the transition criteria based on distance from the 
stagnation point. Thus far, this was not implemented into the Navier Stokes interface in order to avoid 
having to determine a single stagnation point for complex flows. Separating the heat transfer correlations 
from the boundary layer routine allows the Nusselt number to be determined as a point calculation rather 
than integrating from a starting location. This change implicitly assumes that the flow can transition back 
to laminar if the roughness Reynolds number drops below 600 at any location. This change also 
modularizes each step along the process so that developers can change each equation independently. This 
will allow alternate heat transfer correlations to be easily implemented if improvements to this model are 
determined experimentally. 

In addition, GlennICE can use momentum thickness rather than skin friction by using the following 
expression for shear stress: 


T = 


0.1682pV'i 


^ f O 

log 864 -4 ^ 

V V x k 


w 2 


+ 2.568 


JJ 


2t 

pV 2 


Collection Efficiency/Water Mass Flux 

GlennICE uses the same routines as LEWICE and LEWICE3D for determining collection efficiency 
by integrating individual particle trajectories. In addition, the LEWICE collection efficiency module was 
converted to Fortran90 and modularized. This will allow other routines to more easily access the SLD 
models in LEWICE for splashing, bouncing, breakup and reimpingement. Those routines were modified 
such that they are useable in either 2D or 3D. In their current form, developers can also change the 
individual SLD effects by only changing a single routine. In addition, the interpolation routines for air 
velocity were modified to accommodate unstructured grids. A higher accuracy interpolation scheme was 
implemented that matches derivatives at grid boundaries (ref. 55). GlennICE uses the following models in 
the SLD regime: 
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Breakup 


y2^ 

If We d = d — > 13 , then the drop instantaneously breaks up into smaller particles that have the 


same velocity and direction as the original drop. The new drop size is given by 


dyiew — 6.2 d, 


old 


f-f 

K Pa y 


- 1/ 
Re 22 


Splashing/Bouncing 


).859-y/i 


0.859 \Oh w Re, 


1.25 f Pw 


\ 0.125 


If k l = 


V LWC 


(sin 0,j m p ) 


1.25 


> 200, then the fraction mass loss is determined by the 


following equation 


= 0.7 ( 1 - sin a imp ) ( 1 - exp (-0.0092026 ( K L - 200 ))) 

™imp 


The sizes of the splashed particles are given by d sp i ash = 8.72 d imp e 


~° muK where K = Oh w Re|; 25 , 


The velocity and direction of the splashed particles are given by V xsp i ash = V ximp (l.075 - 0.0025 (X ())l/ , ) , 


K 


y, splash 


Vyjmp (o.3 0.002oc im p j and oc sp i as h lan 




A 


y, splash 
\^ x , splash ) 


. A single splashed particle is then 


tracked using the same Lagrangian particle tracking schemes used to determine the initial impact. If the 
particle reimpinges, the splashed mass is added to the new location. 


Ice Density 

For any cell that is rime ( Nf r = 1), the ice density in SLD is given by p/ = 917 kg/m 3 for MVD < 40 pm 
p, = 917-1.5 {MVD - 40) for 40 pm < MVD < 3 1 8 pm and p, = 500 kg/m 3 for MVD > 3 1 8 pm. The 

splashed mass loss and density changes were both determined empirically from previous data (refs. 49 
and 84 to 87). For GlennICE, the ice density correlation is used for both rime and glaze ice shapes. Since 
the correlation defaults to a constant value for MVD < 40 pm, it is used for both SLD = 1 and SLD = 0 
options in the LEWICE and GlennICE input files. 


Energy Balance 

The energy balance equations were re-derived and implemented from scratch in modular form for 
GlennICE. A review of other sources (refs. 56 to 63) determined that there were two differences in the 
formulations used by LEWICE and by other icing codes. The first difference was in the evaporation 
model. While the mass lost by evaporation is very small for an unheated analysis, the energy term can be 
quite substantial since the heat of vaporization is much larger than the heat of fusion. All of the 
formulations use the Chilton-Colbum approximation to tie the mass transfer coefficient to the heat 
transfer coefficient. However there is a slight difference in the implementation of the concentration 
difference. Messinger uses the following equation 
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P V (T S ) P„P r (T,) 


tfevap ^v k m 


T s T a P 

1 P 0 pAL) 

0.622 T 0 T s 


This equation has also been used for scaling conditions (ref. 64) for icing tunnel tests and by Gent 
(ref. 57) in the TRAJICE code after accounting for compressibility effects. LEWICE uses a simplified 
form of this equation by noting that the second term in the denominator is very small for unheated 
surfaces, giving 


k(t.) p„p,(T t ) 


Qevap ' 


1 


0.622 T n 


Myers (ref. 60) and Bourgault (ref. 62) linearize the above equation for easier solution by matrix 
solvers even though the vapor pressure term is highly non-linear. The approximation was justified by 
noting that the surface temperature is usually near freezing even for rime conditions. The vapor pressure 
in LEWICE was derived from a curve-fit of tabulated data and is given as 


ln(P v ) = 20.15247167- 


11097.16963 
1.8 T 


An alternate expression is given by the Goff-Gratch equation (ref. 65) for supercooled water 


log 10 (P v ) = -7.90298 f - 1 


373.16 


l T 


+ 5.02808 logjQ 


373.16 


— 1 .3816 * 10 — 7 1 10 


) 


, 11 . 344 ( 1 - 7 / 373 . 16 ) _ 


+8.1328*10~ 3 |10 


,- 3 . 49149 ( 373 . 16 / 7 - 1 ) 


1 +log 10 (1013.246) 


The Goff-Gratch equation gives slightly higher values at low temperatures. Each of these methods has 
been implemented into GlennICE for evaluation. GlennICE comparisons shown in this report will use the 
formulation in LEWICE. 

The second term that differs in the energy balance is the treatment of the heat lost into an unheated 
surface due to conduction. Messinger ignored this effect. LEWICE assumes an analytical formulation 
based on treating the airfoil as a semi-infinite surface. This gives the conduction loss as 

2k a( T rec- T mp)(Jh-^Jh) 

tfcond ~ I 

y]na a t 2 -t x 

LEWICE can also directly calculate the heat loss due to conduction by using the thermal deicer 
option and turning off all the heat generation. This is an acceptable option to use when interfacing with a 
Navier-Stokes solver, but adds significant computational time otherwise. The direct calculation of heat 
conduction losses is also performed by Myers in the ICECREMO code and by Bourgault in FENSAP- 
ICE. Myers also provides an approximate expression by assuming that no heat is actually lost into the 
airfoil, resulting in a pseudo-steady state heat flux 
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where B is the ice height, which changes with time and T s is the airfoil surface temperature, which is 
assumed not to change. A third approximation can be obtained by solving the heat conduction equations 
analytically using more appropriate boundary conditions. In this formulation, it is assumed that the inner 
surface of the airfoil has an adiabatic boundary condition and that conduction occurs normal to the 
surface. Using a standard separation of variables solution arrives at the following infinite series 


T mp T = ^ 2e ~ a ‘> c " 2t 


T -T 

mp * rec n =Q 


COS 




r c y- 

n y 


with the eigenvalues being given by c n tan(c^) = — 


k:Y 


k a B 


The heat flux at the surface is then 


2k (t —T \ oo 

^ ay rec mp ^ -a c l t - 2 1 \ 

Vcond = “ — L e a “ Cn Sin 2 (c n ) 

1 n = 0 


The first three terms of the series, as well as the Myers conduction expression, have been 
incorporated as an alternate expression for the conduction heat loss. 


Mass Balance 

Myers (ref. 61) provides the following expression for the surface mass balance and film flow: 


dh dB 

Pw-^ + PwV-Q = -Pi-^--PwW- n 


In this expression, W represents the net velocity of water impinging (impingement - evaporation) and 
Q represents the combined effects of pressure, gravity and shear stress. The terms in this expression have 
the following physical meaning 

dh 

— = rate of accumulation 


V • Q = net rate of water transport 
dB 

ft — = rate of ice accumulation 

1 1 


-p w W - n = rate of impingement - rate of evaporation 

However, using the nomenclature used in LEWICE for these terms, this expression becomes 
rh acc^ r hrb= rh i- rh fr- rh e 
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which is the equivalent expression used in LEWICE for the mass balance. The difference between the two 
formulations resides in the value of film height, not in the actual mass balance. LEWICE also uses the 
approximation that film accumulation and runback do not occur simultaneously, hence 


and 


ftv 


dh 

dt 


"ft 


dB 

dt' 


■ -n for t < t s 


PwV • Q = -ft — — Pm ,W n for t > ts 
ot 

where t s is the time needed to reach steady state. The time to reach steady state is determined from the 
approximation for steady film height, hence 


^ss 


PmaxV 


where LEWICE uses = 


a 

ftv I 


W{ 


dV 

ds 


stag 


^ J 


and 

Myers gives h ss — 5 * 10 -6 + 0.40655 ' in 

\ f\T 

The mass balance attributed to Messinger ignored the film accumulation effect and used only the 
steady-state approximation. The mass balance in GlennICE has been formulated such that film 
accumulation and film transport (runback) can occur simultaneously although the Myers formulation has 
not yet been added as an option. There is an additional change that needs to be made to the process for 
GlennICE. In LEWICE, since the ice accretion is 2D (LEWICE3D uses a 2D strip approach for ice 
accretion as well) an explicit marching scheme can be used as the solution mechanism for water film. The 
solution starts at the stagnation point and marches outward from that point. This solution mechanism is 
still valid for 2D flows, provided there is no water film transport into a region where the airflow is 
separated. However, in order to accommodate a full 3D ice accretion process, an alternative solution 
mechanism is needed. The current scheme uses Gauss- Seidel iteration to solve for the water flux. In 
addition, this solution mechanism can also be used for 2D or 2D strip solutions. While the iterative 
solution is slower than the explicit marching scheme in 2D, they provide identical results. The modular 
form of the current scheme was designed to allow for coupling with alternate water flow schemes 
(refs. 66 and 67). Each of these alternative approaches have been added to GlennICE. Due to the large 
number of cases that need to be ran for validation, only the current LEWICE model has been assessed for 
this report. 


Ice Growth 

Once the amount of ice to be added has been determined, the ice volume (or area in 2D) can be added 
to the existing geometry. LEWICE has the capability of growing ice normal to the surface, normal to the 
impacting particle and normal to the flow angle (angle of attack). The latter two options are internal to 
LEWICE and cannot be user-selected. In addition, LEWICE adds or removes points to the surface as 
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needed in order to produce smoother ice shapes. GlennICE can use the LEWICE growth methodology for 
2D geometries, however 3D growth is much more complex. Currently, 3D ice growth can only be 
performed using the three growth direction models mentioned above. The scheme for adding or removing 
points to the surface can only be used in 2D. In addition, Smagglce can only be used on 2D geometries 
and grids must be read in from external files. For these reasons, ice accretion in 3D can only be performed 
in a single time step. 


Experimental Data 

The experimental data described in this paper are the result of a wide variety of tests performed in the 
NASA Glenn Icing Research Tunnel (IRT) in recent years (refs. 38 and 68 to 87). An emphasis was 
placed on recent additions to the SLD database that had not been compared with LEWICE. Eight airfoils 
were selected for this comparison as shown in figure 1. These airfoils and the accompanying ice shapes 
represent the complete set of publicly available data that has been generated in the IRT and digitized for 
single element airfoils. There is some data available on multi-element airfoils, but it was considered to be 
an insufficient amount for validation purposes. There are a total of 1200 IRT runs analyzed for this 
validation report, of which 507 are repeats of previous runs in the IRT. There are 1326 digitized tracings 
at off-centerline locations for a total of 3033 experimental ice shapes. The seven airfoils are as follows: a 
modified NACA 23014, a Large Transport Horizontal Stabilizer (LTHS), a GLC 305 business jet airfoil, 
a second business jet wing, a NACA 0012, a modified NACA 4415, a NACA 0015, and a NLF 0414 



6 

4 

c 2 
>. 0 
-2 
-4 


Large transport horizontal stabilizer 


10 


20 


30 


40 


50 


60 


70 



x, in. 


Figure 1. — Airfoils used in SLD validation. 
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laminar flow general aviation wing. Many of these models were described in an earlier report on the 
LEWICE 2.0 validation effort. In addition to this data set, a newer database has been developed 
consisting almost exclusively of SLD ice shapes (ref. 38). 

The data is taken in the IRT by cutting out a small section of the ice growth down to the airfoil and 
tracing the contour of the ice shape onto a cardboard template with a pencil. The pencil tracing is then 
transformed into digital coordinates with a hand-held digitizer. A flatbed scanner with digitizing software 
has been available to accelerate the data acquisition process. For any given IRT test run, up to five 
spanwise sections of the ice shape are traced and digitized in this manner. There are several steps within 
this process that can potentially cause experimental error. Those that can be quantified by the current 
technique are the spanwise variability, the repeatability error, and errors involved in the tracing technique. 


Results and Comparison Methodology 

This section describes the methodology used to make the quantitative measurements on experimental 
and predicted ice shapes. This methodology has been incorporated into a software utility called THICK 
that calculates and outputs the parameters described. This software was created in order to process the 
large number of ice shapes presented in this report. This program reads two geometry files: one for the 
clean airfoil and one containing an ice shape. This utility has been documented more thoroughly in the 
LEWICE 3.0 User Manual (ref. 29). This software has recently been revised to more accurately capture 
ice features of interest. An additional output file called “ReducedPeaks.dat” provides the location of other 
ice features that may or may not be ice horns that can be evaluated by the user. This file replaces the 
output “Peaks.dat” which often had too many peaks for a user to manually consider when applied to 
experimental ice shapes. The ice shapes are categorized using icing limit, area, horn height, horn 
thickness and leading edge thickness. The parameters were defined in the previous validation report, 
although area is now defined as the true area rather than the integrated thickness. 

The previous section which describes the GlennICE software centered on alternatives to the LEWICE 
model for several features. A comprehensive validation of the software would require validation of each 
feature and combination of features. Time constraints prohibited a comprehensive comparison of this 
type. Instead, validation for the moment will be limited to the GlennICE routines that model the 
same physics used in LEWICE. This limitation results in an exact match for all ice shapes where 
MVD < 40 jam. The modified ice density correlation becomes the only routine that differentiated the two 
programs. This level of agreement between the two software programs was intentional. In order to show 
that a change in a physical model has an effect on ice shape, it was necessary to first have the ability to 
recreate LEWICE results in the new framework. Because of this, results for MVD < 40 jam will not be 
shown as those results were reported on previously (ref. 34). In addition, whereas the discussion below 
may attribute results to LEWICE or GlennICE, the outputs are the same once the density model to 
LEWICE has been modified. Therefore these results are more accurately attributable to LEWICE 3.3 or 
GlennICE 0.1. There are 177 SLD cases to consider in the new database of which 12 are sequential spray 
cases designed to simulate bimodal drop distributions typical of freezing drizzle conditions. In addition, 
there are 521 SLD cases on NACA0012 airfoils from the previous validation of which 132 are sequential 
sprays. 

The comparison methodology will follow previous LEWICE validation efforts (refs. 39 and 40). 
Measured values are the lower and upper icing limit, area, leading edge minimum thickness, lower and 
upper horn thicknesses and lower and upper ice horns. The report on the SLD database (ref. 38) provided 
results for ice height and horn angles. The other parameters (icing limit, area, leading edge minimum) 
were determined by using the THICK utility. The comparison for icing limit is shown in figure 2. This 
figure shows a significant decrease in accuracy for predicting icing limit for SLD cases. While some of 
this difference can be attributed to the difficulty in measuring icing limit experimentally, this factor alone 
cannot fully explain the difference. Much of the difference on the lower icing limit can be attributed to the 
use of a monodispersed drop size as shown by the decrease in the error for the lower icing limit. This 
occurs because the monodisperse MVD results underpredict the icing limit for 75 percent of the SLD 
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Figure 3. — Experimental ice shape with extended icing limit. 


cases. However, use of an 1 1-bin drop size distribution for these cases resulted in a significant over 
prediction for the icing limit. Therefore, it may be necessary to combine a bouncing model that reduces 
impingement limit with the drop size distribution cases to more accurately predict the icing limits. The 
current model just reduces the mass in the downstream region but does not change the icing limit. The 
upper icing limit had an equal number of cases under predicted and over predicted. The differences 
between the code results and experiment are not simply limited to the measured difference in icing limit. 
The experimental ice shapes often have large feathers downstream of the leading edge whereas the code 
results tend to have icing limits that taper off to zero. Figure 3 shows an example of this type of shape. 
LEWICE (and GlennICE) predict the ice shape well past the start of the feather formation, but the three- 
dimensional nature of the SLD feathers is not modeled in either code. Figures 4 and 5 show photographs 
of the lower and upper surface respectively. These figures show that the ice feathers are separated and do 
not form a solid ridge, which is the effect imagined when looking at the two-dimensional tracing. 

Figure 6 shows the comparison for the Leading Edge Minimum Thickness and Ice Area. The ice area 
is actually much better predicted in the SLD regime than for standard drop sizes. This agreement is rather 
fortuitous however. As shown in figure 4, the experimental ice shape in the feather region is largely three- 
dimensional. The tracings, however, are two-dimensional and individual feathers are often not traced. A 
standard experimental technique is used where the pencil is traced over the top of individual feathers 
where it cannot reach to the bottom of each individual feather. This results in the traced experimental ice 
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Figure 4. — Photograph of lower surface ice for case EG1328. 



Figure 5. — Photograph of upper surface ice for case EG1328. 
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Figure 6. — Ice area and LE minimum. 
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Figure 8. — Lower and upper horn angle. 


shape appearing larger than it actually is in the photograph. It can be safely concluded that the codes over 
predict the ice area but the tracing process covers up this result. 

Figure 7 shows the comparison for the lower and upper ice horns while figure 8 shows the 
comparison for lower and upper horn angle. The comparison of SLD shapes to experiment shows 
approximately the same difference as for standard ice shapes. For the experimental ice shapes, the upper 
and lower horns were selected manually rather than allowing THICK to automatically select the largest 
feature. This process was used due to the high variability of the experimental ice shapes. THICK often 
has trouble locating the correct horn locations without manual intervention. Running THICK 
automatically on each shape produced the code values for horn thickness and horn angle. This was done 
to save time in the analysis and was justified on the basis that the code results are smoother and an 
automated process could be used. However, it is possible that THICK did not select the correct horn 
locations for all cases. Figures 9 and 10 show examples of cases with large variations in the upper horn 
location. This result occurs because upper surface feather features were selected for the experimental 
“horn” location whereas the code predicts only the main ice shape and does not model the feather ice 
process. This resulted in a 128° difference in the tabulated measurement for the lower horn angle in the 
first case and a 1 17° difference in the upper horn angle in the second case. Neither result is reflective of 
the general agreement between the experimental ice shape and the code results. Figures 1 1 and 12 show 
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Figure 9. — Example ice shape with large horn angle difference. 



Figure 10. — Example ice shape with large horn angle difference. 


the photographs of the ice shape for these two shapes. In both cases, there are several prominent ice 
features that could be selected as the horn location other than the one selected. 

Figure 13 shows the ice shape that has the best comparison by the quantitative comparison 
methodology. The ice shape comes from a NACA23012 airfoil with the following conditions: c = 12 in., 
V= 200 kts, a = 2°, T 0 = -10 °F, MVD = 225 pm, LWC = 0.55 g/m 3 , t =10 min. This case shows that even 
for a very high water loading, the codes can produce ice shapes that are representative of the experimental 
data. The average comparison is shown in figure 14. The conditions for this case are c = 36 in., V = 

250 kts, a = 2°, T 0 = - 10 °F, MVD = 156 jam, LWC = 0.30 g/m 3 , t = 10 min. The horns on the 
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Figure 11. — Photograph of lower surface ice for case JF1539. 



Figure 12. — Photograph of upper surface ice for 
case I FI 285. 
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Figure 13. — Best comparison from measured parameters. 
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Figure 14. — Average comparison from measured parameters. 
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Figure 15. — Worst comparison based on measured parameters. 
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Figure 16. — Worst comparison No. 1 by qualitative observation. 


experimental shape are slightly longer and closer together than the code results. It is not clear that in this 
case the difference would be aerodynamically significant. Figure 15 shows the ice shape that produces the 
“worst” comparison by the quantitative measures. However, as can be seen from the figure, this ice shape 
is not significantly worse than the previous comparisons. The parameters are worse primarily because of a 
difference in the assigned location for the lower ice horn. The lower horn from the experiment was judged 
to be at v = 4.5 in. while the “lower” horn generated by THICK for the code results was actually on the 
upper surface at x = -0.2 in., y = 0.75 in. While each individual measurement can be defended, the 
outcome results in a measured “difference” of 238° in the reported horn angle, which is not representative 
of the general agreement between the two shapes. It should also be noted that the experimental ice shape 
clearly continues past the end of the template on the lower surface. This causes an erroneous difference in 
the reported icing limit as well. A more valid evaluation of which shapes compare well can unfortunately 
only be performed qualitatively at this time. A utility program was written that creates a Microsoft Word 
(Microsoft Corporation) document containing thumbnail pictures of each shape so that a quick qualitative 
comparison could be made. Figures 16 and 17 show two ice shapes that were judged to be qualitatively 
the worst. The first shape, from a NLF-0414 airfoil where c = 36 in., V= 100 kts, a = 2°, T 0 = -5 °F, 
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Figure 17. — Worst comparison No. 2 by qualitative observation. 

MVD = 45 jum, LWC = 0.97 g/m 3 , t = 10 min, shows that the code can under predict ice feathers on the 
upper surface to a large enough degree such that it undoubtedly impacts the aerodynamics of the two 
shapes. The code was run with a monodispersed MVD for this condition as well. The second shape, from 
a GLC-305 airfoil where c = 36 in., V= 175 kts, a = 6°, T a = -5 °F, MVD = 97 jim, LWC = 0.48 g/m 3 , 
t= 10 min, shows a typical result where potential flow predicts lower ice horns at high angles of attack. 
This type of result is sometimes seen in appendix C cases as well. 


Conclusions 

A new code (GlennICE) has been created for predicting ice accretion on 2D and 3D geometries. The 
intent of this process has been to streamline software development by creating a team-based rather than an 
individual-based development approach which has resulted in a more modularized program in which new 
algorithms can be easily added. The code keeps many of the same features as LEWICE and produces 
similar results for the cases compared thus far. A comparison with 2D experimental ice shape tracings 
confirms the fidelity of the new model. Comparisons in the SLD regime against a new set of experimental 
data confirm the accuracy of the SLD predictions for both LEWICE and GlennICE. A quantitative 
comparison with experimental ice shapes shows that the codes are as accurate in the SLD regime as they 
are in appendix C conditions for horn height, area and horn angle. Agreement is less accurate for icing 
limit. In general, the codes under predict the lower icing limit and show a higher variability in predicting 
the upper icing limit as compared to the agreement seen in appendix C conditions. 

The process used to quantitatively compare ice shape features between two different profiles needs 
further improvement. The results from this and previous comparison efforts have shown that the current 
method can produce misleading results. While visual comparison of ice shape profiles can be used to 
qualitatively assess the degree of comparison, it is difficult to use such an approach to objectively 
determine whether changes in computational simulation methods have resulted in improvement. Currently 
a combination of visual inspection and ice shape feature measurement is the best approach to use for 
validation of computational ice growth simulation methods. 
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